    
    
    rm(list=ls());  options(error=stop)
    source("1.0 MyFunc.R");    source("1.0 Delta_est.R")
    source("1.2 Reg.R"); source("1.3 IPW.R"); source("1.4 G-est.R")
    source("1.5 TR.R")
    source("2.1.1 TSLS.R")
    library(brm)
    MESSAGE <-- TRUE     # compress message (This is a global variable!)
    
    b.index = 1
    
    load("NLS.RData")
    # > colnames(X.mat)
    # [1] "intercept"   "age"         "fatheduc"    "motheduc"    "momdad14"   
    # [6] "sinmom14"    "step14"      "south66"     "black"       "smsa"       
    # [11] "south"       "smsa66"      "kww"         "iq"          "married"    
    # [16] "exper"       "expersq"     "region66"    "fatheduc.na" "motheduc.na"
    # [21] "kww.na"      "iq.na"      
    
    btstrap.each = 500
    
    ### Data preparation
    
    X = X.mat[,c(1:4, 8,9,12, 14, 19,20,22)]
    method   = c("reg","ipw","b-ipw","g","tr","b-tr")
    methods = c("naive", "RDreg", "2sls", method)
    wt = data$weight / sum(data$weight);    wt = wt * N
     
    y.cutoff = median(data$wage)
    
    results.point  = numeric(length(methods))
    names(results.point) = methods
    
    results  = array(dim=c(length(methods),btstrap.each), 
                     dimnames = list(methods,1:btstrap.each))
    names(dimnames(results)) = c("Method", "Bootstrap number")
    
    
    ### Point estimate
    # 0. Naive
    results.point["naive"] = sum(Y[D==1]*wt[D==1])/sum(wt[D==1]) -
        sum(Y[D==0]*wt[D==0])/sum(wt[D==0])
    # 1. 2SLS
    results.point["2sls"] = TSLS(Z,D,Y,X,wt)
    # 2. Proposed methods
    delta.est.list = DeltaEst(X,X,X,X,X,Z,D,Y,wt)
    results.point[4:length(methods)] = delta.est.list$Delta.est
    # 3. RDreg
    rdm = brm(Y,D,X,X, "RD", weights = wt.btstrap, message = FALSE)
    results.point["RDreg"] = mean(tanh(X.btstrap %*% 
                                             rdm$coefficients[1:ncol(X.btstrap),"point.est"]))
    
    ### Bootstrap estimate
    for(btstrap in 1:btstrap.each){
        cat("Bootstrap sample:", btstrap ,"\n")
        Y = as.numeric(data$wage > y.cutoff)
        set.seed(btstrap + (b.index-1)*btstrap.each)
        index = sample(1:N, N, replace=TRUE)
        Z.btstrap = Z[index];  D.btstrap = D[index]; Y.btstrap = Y[index]
        X.btstrap = X[index,]; wt.btstrap = wt[index]
        # 0. Naive
        results["naive",btstrap] = 
            sum(Y.btstrap[D.btstrap==1]*wt.btstrap[D.btstrap==1])/
            sum(wt.btstrap[D.btstrap==1]) -
            sum(Y.btstrap[D.btstrap==0]*wt.btstrap[D.btstrap==0])/
            sum(wt.btstrap[D.btstrap==0])
        # 1. 2SLS
        results["2sls",btstrap] = TSLS(Z.btstrap,D.btstrap,Y.btstrap,X.btstrap,
                                      wt.btstrap)
        # 2. Proposed methods
        delta.est.list = DeltaEst(X.btstrap,X.btstrap,X.btstrap,X.btstrap,
                                  X.btstrap,Z.btstrap,D.btstrap,Y.btstrap,
                                  wt.btstrap)
        results[4:length(methods),btstrap] = delta.est.list$Delta.est
        
        # 3. RDreg
        rdm = brm(Y.btstrap, D.btstrap, X.btstrap, X.btstrap, "RD",
                  weights = wt.btstrap, message = TRUE)
        results["RDreg",btstrap] = mean(tanh(X.btstrap %*% 
                                                 rdm$coefficients[1:ncol(X.btstrap),"point.est"]))
    }
    
    save(results,results.point,delta.est.list,file="./RDataFiles/NLS.btstrap.RData")
    
    
    
    
    
    
    
    
    
